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Abstract 

Theoretical studies of synchronization are usually based on models of cou- 
pled phase oscillators which, when isolated, have constant angular frequency. 
Stochastic discrete versions of these uniform oscillators have also appeared 
in the literature, with equal transition rates among the states. Here we start 
from the model recently introduced by Wood et al. [Phys. Rev. Lett. 96, 
145701 (2006)], which has a collectively synchronized phase, and parametri- 
cally modify the phase-coupled oscillators to render them (stochastically) 
nonuniform. We show that, depending on the nonuniformity parameter 
< a < 1, a mean field analysis predicts the occurrence of several phase 
transitions. In particular, the phase with collective oscillations is stable for 
the complete graph only for a < a' < 1. At a — 1 the oscillators become 
excitable elements and the system has an absorbing state. In the excitable 
regime, no collective oscillations were found in the model. 
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1. Introduction 



The last decades have witnessed the growth of a rich literature aimed 
at developing theoretical tools for understanding the problem of collective 
oscillatory behavior (often loosely termed "synchronization" ) emerging from 
the interaction of intrinsically oscillating units (l|-0|. This development has 
been prompted by the ubiquity of the phenomenon across different knowledge 
areas, with abundant experimental evidence j^-H]. Neuroscience is one par- 
ticular field where this problem reaches the frontiers of our current theoretical 
understanding: neurons are highly nonlinear, interact in large numbers and 
are often subject to noise jsj]. Under these conditions, it is not surprising that 
many theoretical investigations on collective neuronal phenomena (of which 
global oscillations is just a particular case) have been restricted to numerical 
simulations (recent examples include Refs. (j- 

Analytical solutions of the synchronization problem in the presence of 
noise have recently appeared, but at the cost of drastic simplifications of 
each unit. For instance, each stochastic oscillator in the model introduced 



by Wood et al. [121 . 1 131 ] is described by three states connected by transition 



rates, amounting to a discretization of a phase (already a simplification in 
itself @, 0])- The underlying assumption is that the details in the description 
of each oscillator should become increasingly irrelevant as the number N of 



interacting units diverges. This idea was beautifully illustrated in Refs. [12 



13| . which confirmed earlier predictions (derived from a field-theory-based 
renormalization group analysis 14, 15j]) that a phase transition to a globally 
oscillating state belongs to the XY universality class. 

Here we would like to apply this level of modeling to address a different 
problem. We are concerned with global oscillations emerging from interact- 
ing units which, when isolated, are excitable, i.e. not intrinsically oscillatory. 
In a detailed description (say, a nonlinear ordinary differential equation), an 
excitable unit has a stable fixed point (the quiescent state) from which it 
departs to a long excursion in phase space if sufficiently perturbed. It then 
undergoes a refractory period during which it is insensitive to further pertur- 
bations, before returning to rest. A reduction of each unit to a three-state 
system is straightforward in this case 16| : employing the parlance of epidemi- 
ology, each individual stays susceptible (S) unless it gets infected (I) from 
some other previously infected individual, after which it eventually becomes 
recovered (R) and finally becomes susceptible (S) again at some rate. The 
prototypical experimental example of global oscillations of excitable units 
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comes precisely from epidemiology, which shows periodic incidence of some 



diseases [17J (even though a person in isolation will not be periodically ill). 

Models employing such cyclic three-state excitable units have usually 
been termed SIRS models. These can be further subdivided into two cate- 
gories: those with a discrete-time description (cellular automata) and those 
where time is continuous (frequently referred to as "interacting particle sys- 
tems" 18]). In discrete-time models, global oscillations have been reported 19- 



22 



and understood analytically [23] . In continuous-time models, the situ- 
ation is different: in the standard SIRS model with diffusive coupling 24[, 
global oscillations seem not to be stable, even in the favorable condition of 
global or random coupling 24J-I27 . 



From a formal point of view, the cyclic structure of the excitable units 
prevents the use of an equilibrium description. Moreover, the system has a 
unique absorbing state (all units quiescent), which provides an interesting 
connection with an enormous class of well studied problems related to di- 
rected percolation 28||. So, within the context of continuous-time models, it 
is natural to ask: are there markovian excitable systems with an absorbing 
state (i.e. without external drive) which can sustain global oscillations? 

We previously attempted to answer this question by modifying the in- 
fection rate of the SIRS model, which was rendered exponential (instead of 
linear, as in the standard case) in the density of infected neighbors [29[. The 
intention was to mimic the exponential rates of Wood et al. which, in a 
system of phase-coupled stochastic oscillators, did generate stable collective 
oscillations 12|, |l3j. In our nonlinearly pulse-coupled system of excitable 



units, however, this modification actually failed to generate oscillations, but 



rather turned the phase transition to an active state discontinuous [29 . 

Here we make another attempt, now transforming the model by Wood 
et al. so that their oscillating units become increasingly nonuniform. Dif- 
ferently from our previous attempt 29J, units here remain phase-coupled, 



even in the limit where the nonuniformity transforms the oscillators into ex- 
citable units. Only in this limit does the system have an absorbing state, 
and the question is whether sustained macroscopic oscillations survive this 
microscopic parametric transformation. 

The paper is organized as follows. In section [2] we introduce the model. 
Section [3] brings our main results, which are further discussed with our con- 
cluding remarks in section HI 
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2. Model 



Let each unit x (x — 1, . . . , TV) be represented by a phase (j) x , which can be 
in one of three states: <f> x = (i x — l)2n/3, where i x G {1,2,3} (as illustrated 
in Fig. [T]). When isolated, the transition rates for each unit are 

1 — > 2 at rate g(l — a), 

2 — > 3 at rate g, (1) 

3 — > 1 at rate g . 

Parameter a controls the (average) nonuniformity of the uncoupled oscil- 




lators. For a = we recover the uniform oscillators employed by Wood et 



al. |12l, ll3[. For < a < 1, however, each oscillator tends to spend more 
time in state 1 than in the other states, in what would be a stochastic ver- 
sion of a nonuniform oscillation (such as, for example, that of an over damped 



pendulum under the effects of gravity and a constant applied torque |30[). 
For a = 1, an isolated unit stays in state 1 forever (though it will be able to 
leave this state when coupled to other units, see below). The average angular 
frequency of the uncoupled units is u = 2ng (1 — a) / (3 — a), which vanishes 
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continually at a = 1. This is qualitatively similar to what is observed in the 



bifurcation scenario of type-I neuron models [3JJ. Therefore, at a = 1 units 
become excitable. 

The coupling amon g u nits is essentially the same as the one used by 



Wood et al. in Refs. [12l.ll3l]. except that parameter a controls nonuniformity 



in the transition from state 1 to state 2: 



1 _> 2 at rate g 1>2 ^\ = g { e M*»-*)]/* _ , 

2^3 at rate g 2>3 ( ^) = ge^'^, (2) 



k ' k 



3->l at rate g^ , ^ = ^-^/*, 



where g can be set to unity without loss of generality, a is the coupling 
parameter, Ni is the number of neighbors in state i, and k is the total number 
of neighbors. Thus, for a = the original model by Wood et al. is recovered, 
whereas for a = 1 the system has a single absorbing state (all units in state 1). 

3. Results 

3.1. Mean field analysis 

Let gij be the transition rate from state % to state j, given by (j2J), 
where i and j G {1,2,3}. In a mean-field approximation, this leads to the 
following equations: 

A = g3,l( P 3, P l) P S-9lA P l,P2)Pl, (3) 
Pi = gi,2{ P h P 2) P l-92,3{ P 2, P 3) P 2, (4) 
P 3 = 92,3(P2, P 3) P 2-g3,l{ P 3, P l) P B. (5) 

This also coincides with the equations for a complete graph in the limit 
N — > oo with an inherent assumption that we can replace Ni/N with Pj. Nor- 
malization (P3 = 1 — Pi — P2) reduces these equations to a two-dimensional 
flow in the (Pi,P 2 ) plane. 

We start by analyzing the phase diagram of the mean-field equations fl3])- 
(JSJ) in the phase plane (Pi, P2), which is restricted to the triangle ^ Pi ^ 1, 
^ P2 ^ 1 , and ^ Pi + P2 ^ 1 . As shown in Fig. [2J the case of uniform 
oscillators (where the system has a perfect C3 symmetry) shows two bifurca- 



tions as one increases a on the line a = 0. First |l2l . Il3j . at a — a c — 1.5, a 
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Figure 2: (a) Phase diagram with phases characterized by representative phase portraits. 
The triangles mark the border of the physically acceptable region 0^Pi^l,0^P2^1, 
and ^ Pi + P<i ^ 1. Below (above) the dotted line [pink online] the stable fixed 
point is a spiral (node). The remaining lines represent bifurcations in the mean-field 
equations ©-([5]): Hopf (H) bifurcation (red online), saddle-node (SN) bifurcation (blue 
online) and infinite-period (IP) bifurcation (black). The squares indicate the homoclinic 
(HC) bifurcation line. The black triangle indicates the point (a 1 , a') at which the HC, 
H and SN lines meet (see text for details). The dashed line (green online) on top of 
the panel (at a = 1) marks the parameter space region where there is a single absorbing 
state (but not necessarily a single attractor, see text for details), (b) Zoom of panel (a) 
displaying details near the homoclinic line. Panels (c), (d) and (e) show the phase portraits 
for the parameters marked in the panel (b), see also Fig. [4] SLC = stable limit cycle. There 
are stable limit cycles only in the grey regions. 
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supercritical Hopf bifurcation occurs above which the symmetric fixed point 
(P*,P 2 *) = (1/3, 1/3) loses stability and stable collective oscillations emerge 
[a limit cycle in the phase plane (Pi, P2)]. Increasing a further, the frequency 
of these global oscillations decreases continuously as the size of the limit 



cycle increases [32|, [33j . For large values of a, the shape of the limit cycle 
becomes less circular, approaching the borders of the triangle, and global os- 
cillations become highly anharmonic, with a finite fraction of the oscillators 
collectively spending a long time in each of the three states before "jump- 
ing" to the next state at a much shorter time scale. At a = a c ~ 3.102, 
three saddle-node (SN) bifurcations occur simultaneously at the limit cycle, 
corresponding to an infinite-period transition in which C3 symmetry is bro- 



ken (since three stable attractors are created) |34j. As noted previously [32 



this second transition is somewhat artificial from the perspective of more 
realistic models, for which one expects (and observes) oscillators to lock at a 
coupling-independent frequency. However, it is very interesting from the per- 
spective of nonequilibrium phase transitions of interacting particle systems, 
since it provides a spontaneous breaking of C3 symmetry in the absence of 



any absorbing state [34J (as opposed to models with 3 absorbing states, like 
rock-paper-scissors games (35 39]) . 

For a > 0, C3 symmetry is no longer present. In fact, even for arbitrarily 
small nonzero a, the Hopf bifurcation is followed (as a is increased) by an 
infinite-period transition which occurs due to a single SN bifurcation [as 
opposed to the triple SN for a = 0, see Fig.[^a)]. Above this point, collective 
oscillations disappear because units tend to condensate in state 1, which is 
favored by the model nonuniformity [see Fig. [T| . Increasing a further, another 
SN bifurcation occurs [Fig. ^ a)], in which a value of P 2 * > 1/3 becomes a 
second stable fixed point of the model (i.e., state 2, which comes right after 
state 1, now can also attract the oscillators). Finally, for even larger values 
of a, a third SN bifurcation occurs and P 3 * > 1/3 becomes a third attractor. 
With three attractors, the phase space is similar to that observed for a = 0. 

Let us now address the stronger effects of the model nonuniformity. Con- 
sider first the leftmost portion of Fig. ID^a), i.e. for values of the coupling 
a which are too small to yield sustained collective oscillations. In this 
regime, the only attractor is a stable spiral, which for a = is symmet- 
ric: (P 1 *,P 2 *) = (1/3,1/3). Increasing a, this spiral moves towards larger 
values of P*. Eventually, the imaginary part of its eigenvalues vanish and 
the fixed point becomes a node [dotted line (pink online) in Fig. [2^]. Only 
when a = 1 one reaches P* = 1, which is then (and only then) an absorbing 
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state [Fig. EK upper dashed line (green online)]. 

Note that the infinite-period (IP) bifurcation line ends in the beginning 
of the SN and homoclinic bifurcation lines [Fig. EJ^b)]. As we will see below, 
this happens because there is an intrinsic interdependence among these three 
bifurcations. An infinite-period bifurcation is just a SN bifurcation in which 
the saddle and the node are born at the limit cycle, whereas in the homoclinic 
bifurcation the limit cycle collides with a pre-existing saddle and becomes a 
homoclinic orbit. In both cases the period diverges at the bifurcation, though 



with different scaling regimes [30 



Now we turn to the main question of this article, namely, what happens 
if we start from a regime with collective oscillations (say, with a > a c , a > 0) 
and increase the system nonuniformity? We will see that many scenarios 
exist, but all of them ultimately lead to the destruction of the limit cycle 
for some a ^ a' < 1 [Fig. Et^b)]. We consider first the simplest cases. If, 
on one hand, we fix a value of a very close to a c , the size of limit cycle 
will continuously decrease with increasing a to vanish in a supercritical Hopf 
bifurcation [H line (red online) in Figs.^a) andEl(b)]. If the fixed value of a is 
sufficiently close to a c , on the other hand, the limit cycle will be destroyed by 
increasing a via an infinite-period bifurcation [IP line in Figs. |2]^a) and[2^b)]. 

There are intermediate values of a for which the way to the annihilation 
of the limit cycle is more tortuous [Fig. E^b)]. The key player to watch 
in these scenarios is a point (a', a') (with a c < a' < a c and < a' < 1) in 
parameter space where the Hopf bifurcation line ends [represented by a black 
triangle in the Figs. Eta) andEl(b)]. For any a G (a c , a'), increasing a will lead 
to extinction of the limit cycle via a Hopf bifurcation [Figs. EJ^a) and (b)]. 
For a J$ a', the limit cycle is first destroyed by an infinite-period bifurcation, 
then emerges again through a homoclinic bifurcation, and finally disappears 
in a Hopf bifurcation [Fig. El^b)]. There is also a small reentrance in the 
homoclinic bifurcation curve [Fig. E^b)], so that for any value of a in this 
short interval, the growth of a leads to destruction and recreation of limit 
cycle by two consecutive homoclinic bifurcations. Finally, for smaller values 
of a, a SN bifurcation leads to the phase portrait of the Fig. El(d) before the 
limit cycle is annihilated by a Hopf bifurcation [Fig. El(b)]. 

If we set a > a' and a < a', we will have only one stable node with P* « 1 
(remember that for a = 1 we have P* — 1, which corresponds to the absorb- 
ing state). Increasing a, there will be a first SN bifurcation at which a saddle 
is born with an unstable node (which quickly becomes a spiral) [Fig. Et^a)]- 
Further increases in a will give rise to two additional SN bifurcations, as 
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described above. After the third SN bifurcation, the system will have three 
attractors with a phase portrait of the same type as that of the figure [3(b). 
Importantly, besides the three saddle-node bifurcations mentioned, no other 
bifurcation occurs for a > a'. In particular, we found no limit cycles for 
a' < a < 1, so that the nonuniformity alone can destroy collective oscilla- 
tions, even in the absence of an absorbing state. 

3.2. Complete graph simulations 



As mentioned in subsection I3.1[ the mean-field calculations are exact for 
the complete graph in the thermodynamic limit. For completeness, we now 
study the effects of finite-size fluctuations. While we do not expect stable 
collective oscillations to appear for a ^ a', stochastic oscillations have been 



predicted to appear around (mean-field) stable spirals [40| in finite popula- 
tions. Moreover, fluctuations become particularly prominent in systems with 
multistability, ultimately determining the attractor to which the system will 
converge. 

First, let us deal with the effects of finite-size fluctuations in a scenario 
with multistability by starting the system near an unstable stationary so- 
lution at the intersection of the basins of attraction of the three stable so- 
lutions [Fig. [3(b)]. At point (b) in Fig. [3(a), tristability is revealed by the 
different fates of the system for three independent realizations of the same 
macroscopic initial condition. 

Next, after each sample was already in its respective steady state, we 
have increased a discontinuously [(b)— >{c) in Fig. [3(a)] without interruption 
in the simulations in order to extinguish one of the attractors. Therefore, the 
sample which was around the newly extinct attractor converges to the stable 
solution of the new larger basin of attraction [while the other two samples 
are only slightly disturbed, see Fig. [3(c)]. 

Then, similarly to the previous case, we decrease a discontinuously [(c)— >(d) 
in Fig. [3(a)] in order to suppress the attractor of the isolated sample, which 
is forced to converge to the remaining stable fixed point. At this stage, all 
three samples have the same behavior, apart from fluctuations [Fig. 13(d)] . 

Finally, we reduce a and a discontinuously and simultaneously [(d)— >(e) 
in Fig. [3(a)] in order to eliminate the last stable fixed point. The three 
samples thereafter converge (at about the same time) to a stable limit cycle 
in which they oscillate almost together, apart from fluctuations [Fig. [3(e) 
and (f)]. Note that the process illustrated in Figure [3J does not close a cycle 
and is irreversible: no matter whether we change the control parameters to 
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Figure 3: Phase portraits and simulation results for a complete graph with TV = 500. Every 
50 time units, the parameters of the model are changed so as to eliminate a stable fixed 
point of the model, (a) Phase diagram with labels to indicate the parameter-space places 
of the phase portraits (b)-(e) and arrows to show the order of the parameter-value jumps, 
(b) Phase portrait for a — 0.2 and a = 4. The three samples of the simulations begin 
at the unstable fixed point close to triangle center, then each one converges to a different 
attractor owing to fluctuations, (c) Phase portrait for a = 0.9 and a = 4. During the 
simulations, we change the parameter a to the value of this phase portrait. Thus, one of 
the attractors disappears, so the sample which was around one attractor converges to the 
stable fixed point of the attraction basin where this sample is now. Obviously, the other 
two samples continue around their attractor. (d) Phase portrait for a = 0.9 and a = 2.5. 
Similarly to panel (c), another attractor disappears. Thus, the sample which was isolated 
converges to the remaining stable fixed point, (e) Phase portrait for a = 0.2 and a = 2.25. 
Finally, the remaining stable fixed point disappears, so the three samples oscillate almost 
together (except for fluctuations). Panel (f) shows time series for simulated trajectories 
shown in panels (b), (c), (d), and (e), in chronological order. The false impression that 
the inequalities < P 3 < 1 arc sometimes violated is caused by the thickness of the lines. 



their initial values [i.e. (e)— >(b) in Fig. E^a)], or we apply the reverse changes 
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Figure 4: Phase portraits and simulation results for a complete graph with TV = 20000, 
sec Fig. [2j (a) Phase portrait for a — 0.69 and a = 1.675. There are two attractors: a 
node and a spiral (they are separated by a saddle), (c) Phase portrait for a = 0.69 and 
a = 1.725. The spiral becomes unstable via a Hopf bifurcation and is surrounded by a 
stable limit cycle, (e) Phase portrait for a = 0.69 and a = 1.775. After a homoclinic 
bifurcation, the limit cycle disappears. Panels (b), (d), and (f) show time series for 
simulated trajectories, respectively, shown in panels (a), (c), and (e), with examples of 
collective excitability (dashed, blue online) and stochastic oscillations [arrow in the inset 
of (b)]. 
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[(e)— >-(d)— >(c)— >(b)], it is extremely unlikely that the system will return to 
its initial steady state [Figure 0(b)] . 

Now consider the phase diagram of Figure |2j(b). In the far left (lower 
values of a) there is only one stable node. Increasing a, a SN bifurcation 
occurs in which a newly arisen node turns into a spiral after an extremely 
small increase in a. Thus, we get the phase portrait of the Figure [5Jc), 
whose simulation results for a complete graph appear in Figures |U(a) and (b). 
Although the mean-field calculations predict that the spiral in Figures |2](c) 
and |4](a) is stable, even for large complete graphs (N = 20000 in Figure H|) 
finite-size fluctuations eventually lead to system to the other (fixed point) 
attractor. However, before that happens, these fluctuations lead to stochastic 



oscillations around the spiral [40| [see arrow in the inset of Figure H£b)]. 

Increasing a further, the spiral loses stability in a Hopf bifurcation, giv- 
ing rise to a stable limit cycle (SLC). We arrive then at the phase portrait 
shown in Figure EJ^d), with simulation results for a complete graph shown in 
Figures HJ^c) and (d). Again, while the mean-field calculations predict that 
the limit cycle in Figures E|d) andHJ^c) is stable, finite-size fluctuations lead 
the system to the other attractor in a considerably shorter time than in the 
previous case [Figs. M^c), @Ja) and (b)]. 

Soon after the Hopf bifurcation, increasing a leads to an increase in 
the size of the SLC, until it is finally destroyed by a homoclinic bifurca- 
tion [Fig. [2](b)] . Thus, we reach the phase portrait of the Figure E£e), with 
the results of simulations of a complete graph shown in Figures |U(e) and (f). 
Now there is only one attractor in the system. As expected, if we start the 
system in the unstable spiral, there will be oscillations of increasing ampli- 
tude until the system reaches the stable steady state. 

Note that the presence of the saddle in Figs. |2|c)-(e) sets the conditions 



for collective excitability [29|. If the system is initially in the stable node, 
perturbations that throw the system to a different point still within its basin 
of attraction can lead to qualitatively different relaxation processes. Small 
perturbations decay monotonically, whereas large perturbations will take the 
system further away from the fixed point before returning to it (because the 
system is required to go around the stable manifold of the saddle). This can 
be clearly seen in the simulations (compare the dashed and thick lines - 
blue and pink online, respectively — in Figure H]). 
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4. Concluding remarks 



We have proposed a variant of the model by Wood et al. 12|, [l3| in which 
stochastic oscillators can become increasingly nonuniform as parameter a 
goes from to 1. We have obtained the phase diagram of the mean-field 
version of the model in the (a, a) parameter plane. In addition to the previ- 



ously reported phase transitions for a = [12|, |13|, l34| , we have obtained for 
q^O several bifurcations in the mean-field equations of the model, including 
saddle-node, infinite period, Hopf and homoclinic. Collective excitability 29 
has been shown to occur in some parameter regions, as confirmed by sim- 
ulations of complete graphs. Simulations have also confirmed the overall 
predictions of the mean-field analysis, although the stability of some sta- 
ble limit cycles and fixed points has failed to resist the effects of finite-size 
fluctuations. 

In the regime in which the units are excitable elements, we did not find 
stable global oscillations for the particular choice of a nonlinear coupling de- 
fined by Eqs. |2J even in the complete graph. This topology is the one in which 
a synchronized solution would be most favorable, as shown in a number of 
previous reports [lil 



13, 19, 32, 33, 40 



which retrospectively justifies why 
we have not attempted to run simulations of the model in a hypercubic (or 
even small- world) topology. We have shown that, in the model of Wood et 



al. Il2j,[l3|], nonuniformity hinders the synchronization among the oscillators, 
while large enough nonuniformity and, consequently, the transformation of 
oscillatory units into excitable elements prevents the emergence of a synchro- 
nized stable state. It remains to be investigated whether this can be achieved 
with a different type of coupling. 

Our results nonetheless raise some interesting questions which deserve 
further investigation. For any nonuniformity (a > 0) the transition to a 
synchronized state in the complete graph occurs in the absence of C3 sym- 
metry. Does the transition still occur in hypercubic lattices under these 
conditions 12]? If so, do the critical exponents depend on a? Investigations 
in these directions would certainly be welcome. 
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